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We consider the NMR signal from a permeable medium with a heterogeneous Larmor frequency component 
that varies on a scale comparable to the spin-carrier diffusion length. We focus on the mesoscopic part of the 
transverse relaxation, that occurs due to dispersion of precession phases of spins accumulated during diffusive 
motion. By relating the spectral lineshape to correlation functions of the spatially varying Larmor frequency, we 
demonstrate how the correlation length and the variance of the Larmor frequency distribution can be determined 
from the NMR spectrum. We corroborate our results by numerical simulations, and apply them to quantify 
human blood spectra. 



INTRODUCTION 

Transverse relaxation of the NMR or ESR signal acquires 
distinctive features [1] when it occurs in a medium which pos- 
sesses magnetic structure on a mesoscopic scale. This scale 
is intermediate between the microscopic atomic or molecu- 
lar scale, and the macroscopic sample size or the resolution 
achievable with time-resolved Faraday rotation measurement 
[2], magnetic resonance imaging [3] (MRI), or ESR imaging 
[4]. In this work we consider the transverse relaxation in a 
broad class of media where the mesoscopic structure can be 
characterized by variable Larmor frequency f2(r). 

The problem of mesoscopic contribution to transverse re- 
laxation arises in a broad variety of contexts, ranging from 
solid state physics to radiology. In the semiconductor spin- 
tronics, the spatially dependent Larmor frequency ft(r) for 
electrons or holes can be induced either via ferromagnetic im- 
printing, or electrostatically by locally varying the electron 
g-factor [5-9]. For nuclear or electron spins in liquids, the 
spatially dependent Larmor frequency can arise due to hetero- 
geneous magnetic susceptibility. The latter property is crucial 
in the field of biomedical MRI, where the heterogeneous sus- 
ceptibility x( r ) is inherent to majority of living tissues due 
to paramagnetism of deoxygenated haemoglobin in red blood 
cells [10, 11], enabling in vivo visualization of regional activa- 
tions in human brain via functional MRI [12-15]. Moreover, 
the susceptibility contrast can be enhanced by doping blood 
with magnetic contrast agents for clinical purposes, e.g. for 
diagnostics of acute stroke. 

The common property of all these systems is the presence 
of the transverse relaxation that occurs via the dispersion of 
precession phases of spins accumulated during their diffusive 
motion. This reduces the vector sum of magnetic moments 
from all spins in the sample and causes attenuation (termed 
dephasing) of the measured signal time course s(t). The sig- 
nal s(t) is typically a sum over a large number of spins ac- 
quired over a macroscopic volume V whose size greatly ex- 
ceeds the diffusion length Id oc \ft. The randomness of the 
Brownian trajectories results in an effective averaging (diffu- 
sion narrowing) of the contributions from different parts of 
the volume V that are separated by less than Id- Further av- 



eraging occurs between larger domains of size exceeding Id- 
This self-averaging character of the measurement is a major 
obstacle in quantifying the properties of the medium on the 
scales below V. Here we study how the geometric structural 
details of f2(r) on the mesoscopic scale are reflected in the 
measured signal, i.e. survive the above averaging. For defini- 
tiveness, we will use the established NMR terminology, since 
generalizations onto the ESR case present no problem. 

Theory of mesoscopic relaxation in the presence of the dif- 
fusion narrowing has been previously addressed in the MRI 
context [16-24]. These works, however, either do not di- 
rectly relate relaxation to the structure [17], or employ strong 
simplifying assumptions (representing the medium as a dilute 
suspension of mesoscopic objects with small volume frac- 
tion ( < 1 [16, 18-24] and a particular shape [16, 18- 
21, 23, 24]). This leads to the virial expansion for the sig- 
nal s(t) oc e - ^** 1 , where f(t) describes the dephasing 
due to a single object. The mesoscopic field heterogene- 
ity is accounted for perturbatively, which yields relaxation 
/ oc f2 2 . As a consequence, the known results in the diffusion- 
narrowing regime are limited to dilute suspensions of effec- 
tively weakly magnetized objects. Examples are dilute blood 
samples, or tissues with small volume fraction of paramag- 
netic vessels, in the fields Bq < 1 T. 

The aim of this paper is to provide a framework for the 
transverse relaxation from arbitrary magnetic media beyond 
the limitations of dilute suspension and the weak magnetiza- 
tion. First, we suggest a universal description of the spectral 
lineshape s(oj) of the signal: 



where function is a measurable characteristic of the 

medium. The quantity —Yi(u>) can be loosely interpreted as 
a dispersive relaxation rate, with Eq. (1) being a generaliza- 
tion of the conventional Lorentzian line shape to the case of 
heterogeneous magnetic medium. Technically, the represen- 
tation (1) originates from the standard form of single particle 
Green's function in many-body physics, with S(ijj) called the 
self-energy part [25], the term we will use here. The imme- 
diate advantage of using instead of the traditional line- 
shape s(u>) is that the frequency-dependent part of the self- 
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energy is entirely determined by the mesoscopic magnetic 
structure and thus vanishes for homogeneous media. In con- 
trast, the effect of mesoscopic structure on the lineshape s(to) 
is a complicated deformation of a Lorentzian, that is difficult 
to quantify [26, 27]. Further we relate the self-energy to the 
geometric structure on the mesoscopic scale for media fully 
permeable for spin carriers, and verify the results using ab in- 
tio simulations of the transverse relaxation. Finally, we apply 
our general results to quantify the line shape of water proton 
resonance in blood [27] in terms of mesoscopic structural pa- 
rameters. 



RESULTS 

Here we consider the NMR signal from a random 
medium which is characterized by the correlation functions 
r„(ri, r„) = (0(ri)...0(r„)). These functions vary on 
the mesoscopic scale much smaller than the size of the acqui- 
sition volume V. The self-averaging character of the measure- 
ment implies that one needs to find the signal from a particular 
realization of il(r), and then to average over the realizations 
according to the distribution moments T n . In what follows, 
we suppose for simplicity that the medium is isotropic and 
translation-invariant in a statistical sense. These assumptions 
cover a wide variety of applications, such as the transverse 
relaxation in blood, or in the brain gray matter. 



Spectral lineshape 

A spin traveling along Brownian path r(t) originating at 

r(0) = r , acquires the relative phase exp 

Here O(r) is the variable component of the Larmor frequency; 
(O(r)) = 0. The time evolution 

S(r ,r;t) = f Vr(t) e~ { Jo'^K')]-/o'^ 2 /« (2 ) 

of the magnetization packet initially at r = ro, 
<?(r , r, i)|t = o = S( r — r o)> is the precession phase averaged 
over the Wiener measure on the diffusive paths. Equivalently, 
Q(ro,r;t) is the Green's function (fundamental solution) of 
the Bloch-Torrey equation [28] 



d t tp = V r (r>V r V) -ift(r)V>. 



(3) 



The diffusive dynamics (3) leads to the transverse relaxation, 
^(ro.r,*)!^ =0. 

For simplicity, and in order to underscore the effects of sus- 
ceptibility contrast, below we assume homogeneous diffusiv- 
ity D(y) = D = const. We also excluded conventional expo- 
nential factors associated with homogeneous components of 
Larmor frequency and of transverse relaxation. 

The mesoscopic component s(t) of the measured MR sig- 
nal (normalized to s| i=0 = 1) is obtained in two stages. First, 



one considers the signal sq from a given realization of 17 (r). It 
is given by the magnetization Q(r ,r; t) that is summed over 
the final positions r which spins reach during the time interval 
t, and averaged over initial spin positions r . The signal s(t) 
is then found by averaging of so over all realizations of O(r), 

*(*) = <*>(*)>= 1^ (g(r 0l r;t)) = G(f;k)| k=0 . 

(4) 

Here we used the translation invariance of the distribution- 
averaged Green's function (Q(ro,r;t)} = G(r — ro,t), 
with its Fourier components defined as G^.k = 
JdtdrG(r,t)e i0Jt - ikr . 

Eq. (4) connects the spectral lineshape s(ui) = G Wj k|k=o 
with the distribution-averaged Green's function of Eq. (3). 
The result of such an averaging (see Methods) may be rep- 
resented as 



where the diffusion propagator 
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(5) 



(6) 



-iw + Dk 2 

is the Green's function of Eq. (3) with O = 0. 

The generic representation (5) yields the dynamics u> = 
uj(k) as the pole of the propagator G^,k, modified by inter- 
action with complex environment embodied in S Wi k [25]. For 
the ballistic dynamics, expansion of £ Wj k would give the re- 
fraction index of the medium (when u oc k), or the renor- 
malized quasiparticle mass and lifetime (when u oc k 2 ). In 
the present case of the diffusive dynamics (3) (iu> oc fc 2 ), 
the self-energy S w k contains all the measurable informa- 
tion about mesoscopic relaxation. In particular, its expansion 
S w ,k = S(0) — (5D)k 2 + ... in u> and k describes the effect 
of the mesoscopic relaxation on the coarse-grained dynamics 
of the magnetization density. In other words, Eq. (5), viewed 
as an operator relation, is the effective Bloch-Torrey equation, 
that acquires higher derivatives in t and r, after averaging over 
the medium. According to Eq. (4), the signal lineshape (1) is 
a particular case of (5) with 



S(w) = S w; k|k=0 • 



(7) 



Weak dephasing: A perturbative solution for the lineshape 

As a simplest example we now consider the weakly mag- 
netic medium whose Larmor frequency dispersion Sfl = 
\J (f2 2 ) is small, and find the relaxation in the lowest order 
in (<5fi) 2 = r2(r)| r= o. The self-energy is given by the single 
Feynman graph (see Methods, Fig. 4(c)), 

d d q T 2 (q) 



v^pcrt 
_2 Vk 



(8) 



(2n) d -iu; + D(k + q) 2 ' 

The lowest order contribution to the self-energy part for the 
signal, 

1%) 



V ' J (27T) d 



—ioj + Dq 2 



(9) 
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involves only the angular-averaged two-point correlation 
function ^(g) = (r2(q))q. Eq. (9) is the Gaussian approxi- 
mation in a weakly paramagnetic medium. Corrections to this 
approximation, involving correlators r„, first appear in the n- 
th order in the Born series for £ w ,k- 

Equation (9) gives the universal short-time asymptotic be- 
havior: In the limit ui — > oo this is the leading in 1 /uj contribu- 
tion to £(<x>), yielding ^(ui^^^oo — (Sfl) 2 /iu>. Substituting 
it into Eq. (1) and taking the inverse Fourier transform, results 
in s{t)\ t ^ ~ 1 - (5CI ■ i) 2 /2, for 5Q ■ t < 1. 

Connection to the susceptibility structure: Locality 

The diffusing spins sense the distribution of the Larmor fre- 
quency offset f2(r). Often times, in paramagnetic media, such 
a distribution is induced by a variable magnetic susceptibil- 
ity profile x( r )> which is practically interesting. Then the 
medium is characterized by the correlation functions of x, as- 
suming zero average (x) = 0. In particular, here we consider 
the two-point correlation function (r) = (x(r)x(O)) that 
for a statistically isotropic medium depends only on r = |r|. 
For nonferromagnetic media, such as deoxygenated blood, 
|x| <C 1, hence the induced field Q(r) is connected to x( r ) 
by the convolution with an elementary dipole field Y(r), 

fXr) = 47rfWr)* X (r), Y(r) = ^(^--lj. 

(10) 

Here fl is the uniform Larmor frequency component, and 
Y(r) includes both the local contribution oc <5(r), and the 
Lorentz cavity field, that compensate each other in three di- 
mensions [29, 30]. 

The angular-averaged correlators of Larmor frequency and 
of the underlying susceptibility are proportional to each other 
(locality [22]): 

I^(r) = (n(r)fi(0)) f = (47rAn ) 2 -r^(r), (11) 
where A 2 = JdkY(-k)Y(k) . (12) 

Here Y(k) is the Fourier transform of the dipole field; in d = 
3 dimensions, Y(k) = 1/3- k 2 z /k 2 , and Eq. (12) yields A 2 = 
4/45. In d = 2 dimensions, F(k) = 1/2 - fc 2 /fc 2 , yielding 
A 2 = 1/8. 

To prove Eq. (1 1), we note that the Fourier transform of the 
field ft induced by the susceptibility profile with the orien- 
tation TZ is 0-R.(k) = 47rfi F(k)x(7?.~ :L k). Substituting tt n 
into first Eq. (11) and averaging over orientations TZ yields the 
right-hand side. Physically, we went from averaging over the 
orientations TZ of the medium relative to a given z-direction of 
the main field Bqz to averaging over the direction of the field 
at fixed orientation of the structure x( r )- 

The locality property (1 1) means that, in the second order in 
fl, after the orientational averaging the diffusing spins effec- 
tively interact directly with the susceptibility profile x( r )> Via 




FIG. 1: The two-dimensional medium, (a) The susceptibility profile x( r ) 
generated by the random self-avoiding addition of disks of radius p with 
volume fraction £ = 0.461, and the corresponding Larmor frequency off- 
set Q(r) induced by a vertically applied field. The shown fragment is 16 
times smaller in each dimension than the whole simulation box of size 300p. 
The diffusion of spins was modeled as random hopping on a square lattice 
(4096 X 4096) with lattice constant 300p/4096, in the field shown in right 
panel. Time was measured in units of t p = p 2 /D. (b) Correlation func- 
tions rjj-(r) and T2(r) of the susceptibility and of the Larmor frequency. 
The images are zoomed 4-fold as compared with (a), (c) Illustration of lo- 
cality in d = 2. Left panel: Coinciding angular-averaged correlators of the 
susceptibility r¥(r) (red) rescaled by the factor = 1/8, and of the Lar- 
mor frequency (blue). Dashed horizontal line corresponds to expected value 
f(l — C)/8 atr = 0. Right panel shows the angular-averaged Fourier trans- 
form of the correlator, T%(k), with the pronounced peak at k c fa 2.3/p. 
Noise increases for small k due to finite size effects; T^lfc— o = 0. The 
mean k of the first peak found with the weight fe, which is inherent to the 2d 
integration, is k c = 2.20/ p. 



47rAS!ox(r). This interaction is much simpler than that with 
the susceptibility-induced field fi(r), which involves a convo- 
lution with the non-local field (10) with a complicated angular 
dependence. The underlying reason for this fortunate property 
is the scaling Y(r) oc r~ d in d dimensions [22]. The local- 
ity (1 1) is demonstrated for the two-dimensional numerically 
generated medium in Fig. 1. 
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Model of magnetic structure 

Aimed at relating the signal to the magnetic structure, we 
consider the medium, embodied by the distribution of x( r )> 
that has a single mesoscopic length scale l c . This scale, and 
the dispersion SSI, are the two parameters that we employ here 
to characterize the medium. 

For the model calculations below we focus on the angular- 
averaged correlation function r 2 entering Eq. (9). The pres- 
ence of a single length scale leads to a well-defined dominant 
peak in the function r 2 (fc). Such a form is demonstrated for 
the d — 2 case in Fig. 1(c). We suggest to approximate such a 
peak via a delta-function 

d = 3: r^(jfc) = \{8Vt) 2 l 2 c 5{k-k c ) , k c = 2n/l c , (13) 

normalized to / ^J^ d T 2 = (Sfl) 2 , assuming that the contri- 
bution of other length scales is less relevant. In real space, the 
correlator (13) 



sin k c r 
k c r 



(14) 



In d = 2 dimensions, the corresponding correlator is r 2 (k) = 
(SQ) 2 l c S(k - fc c ) and T^(r) = (<5ft) 2 J (k c r), where J 
is the Bessel function. Generalization onto the case when 
the medium has a set of harmonics (14) with different k c is 
straightforward: e.g. for d = 3, 

T 2 ~{k) = ) 1 {5nfY,Pjlc)S(k-k cj ) with $>, = 1. 

3 3 

(15) 

With the approximation (13), evaluation of self-energy (9) 
in the lowest order is trivial: 



-£ pert H 



a 2 /t c 



1 — iut c 

Here we introduced the diffusion time 

t c = 1/Dk 2 = (l c /2ir) 2 /D 



(16) 



(17) 



past the Larmor frequency correlation length l c = 2ir/k c . 
Time t c sets the scale for both the frequency cj and for the 
self-energy. The dimensionless perturbation series parameter 

a = SQ-t c (18) 

controls the dephasing strength. 

The lineshape (1), (16) is a sum of the two Lorentzians 



Wi W2 



s(uj) -- 

ei — iuj e 2 — iuj 

s(t) = Wl e- eit + w 2 e~ f - 2t 



Wl +w 2 = l, (19) 



where (keeping terms up to the order a 2 ) 



eit c = 1 — a 2 , e 2 t c = a 2 . 



(20) 



(21) 



and the weights w\ = —a 2 and W2 = 1 - 
(9) is valid for a <C 1. The signal s(t) 



a . The approach 



1 



W{t/t c f 



for t t c , and decays monoexponentially for t — > 00 as 
s ~ w 2 exp{-t/T 2 pcrt } with the rate 1/T 2 pcrt = a 2 /t c = 
(5£l) 2 t c . The latter rate arises due to angular diffusion of 
the precession phase: As the phase changes by <~ a at 
each time step t c , over large time t ^> t c the rms phase 
~ ayjtjtc reaches ~ 1 fort ~ j-P crt (Dyakonov-Perel relax- 
ation [31]). The relaxation rate decreases with decreasing l c , 
as the medium effectively becomes more homogeneous (dif- 
fusion narrowing). 



Strong structural dependence 

What happens in a complex medium with a few harmon- 
ics, Eq. (15)? In the lowest order, O (((5£!) 2 ), harmonics with 



t C j con- 



different l C j and dephasing strengths at- 
tribute to the self-energy (9) additively. Thus S(w) is a sum of 
n terms (16), and the lineshape is a sum of n + 1 Lorentzians 
that correspond to the presence of n distinct length scales. 

We emphasize that the apparent n + 1-exponential form of 
the perturbative result s(t) (biexponential (20) for n = 1) here 
has nothing to do with existence of n + 1 macroscopic "com- 
partments" with different relaxation rates. (Note also that the 
weight w\ < in Eq. (20).) This simple example is a warning 
against common practice of literally interpreting bi- or multi- 
exponential fits. 

The t — > 00 relaxation is determined by the pole of sig- 
nal (1) closest to real axis. For weak dephasing, such a pole 
is given by i£ pcrt (0), yielding additive contributions for the 
total rate 1/T 2 pcrt = £)™=i a 2 /t cj similar to the Matthiessen 
rule in kinetic theory. Equivalently, 



1 



ypcrt 



/ 



d 3 rdV (n(r)n(r')) 
V 47rL>|r-r'| 



(22) 



(This result can be already deduced from Eq. (9) by setting 
there u = 0.) 

The rate (22) is formally equivalent to the Coulomb en- 
ergy of a fictitious charge distribution with the charge den- 
sity O(r) [or of 47rAf2ox(r) via locality]. Extending this 
mapping, the variable diffusivity D(r) would be analogous 
to a variable dielectric constant. Equivalence of (22) with 
Coulomb problem is due to the Laplace form of Eq. (3). In 
the limit t — > 00, relaxation is determined by the diffusing 
spins wandering infinitely far to explore the magnetic struc- 
ture. The spins then become mediators of the effectively long- 
ranged interaction between the different parts of the magnetic 
structure. The long-time asymptote corresponds to the time- 
averaged diffusion propagator which becomes a Coulomb po- 
tential, f^dt G°(r, t) ocl/r. 

The self-energy (9), as well as the rate (22), strongly de- 
pend on the geometric structure, as the convergence of inte- 
grals is determined by the specific way of how T 2 vanishes 
at short distances (large k). This nonuniversality can be also 
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seen from mapping onto Coulomb problem, since XlP ort maps 
onto capacitance which is sensitive to the conductor geometry. 



Result for lns(i) 

The self-energy (9) substituted into the lineshape (1) is a 
natural extension of the results [22] (rederived in a different 
way in [23, 24]) to the case of arbitrary volume fraction (. 
The signal takes the form s(t) = exp[— /(f)], where 



whose solution is 



/H - -- 



f d d q 



r 2 (q) 



{w + iO) 2 J (2Tr) d -ioj + Dq 2 



(23) 



Indeed, while f(t) <C 1, Eq. (23) is the lowest order expan- 
sion of (1) in E, valid for t « t c /a 2 . For t > t c , Eq. (23) 
yields the same relaxation rate as Eqs. (1) and (9). As these 
two time domains overlap, the equivalence is proven for all t. 
Relaxation in the dilute suspension is just a particular case of 
(23), with T 2 ~ Cr^. where is the correlator for a single 
object. 



Moderate dephasing: The self-consistent Born approximation 
for the lineshape 

Below we attempt to go beyond the perturbative calcula- 
tion and explore the case of moderate dephasing, a < 1. An 
exact solution would amount to summing up all the diagrams 
for the self-energy to all orders in a, an arduous task even 
within a simplifying assumption (13). Here we consider the 
self-energy in the self-consistent Born approximation (SCBA) 



_ y = f d d <i I^q) 

J {2v) d -iuj + D(k + q) 2 - E w>k+q • 



(24) 



Technically, Eq. (24) amounts to summing all the non- 
crossing Feynman graphs (see Methods, Fig. 4). The SCBA 
is akin to mean field theory (more precisely, it is equivalent to 
the mean field treatment of the nonlinear term that arises after 
Gaussian disorder-averaging in the replica field theory). Its 
advantage is that one can move quite far analytically, by sum- 
ming up an infinite subset of Feynman graphs, thereby collect- 
ing contributions from all orders in the external random field 
O(r). Its disadvantage is that it is uncontrolled since it leaves 
out an infinite subset of graphs for S w k . 

Eq. (24) is a complicated integral equation for E w ,k due 
to the k-dependence of S on the right-hand side. Since we 
really only need E|fc = o, we now make another (uncontrolled) 
simplification that will lead to an ansatz for E(w). Neglecting 
the fc-dependence of the self-energy in the denominator, we 
obtain the quadratic equation 



-E( W ) c 



E( W ) = 



1 — iujt c — R(to) 
2E 



-IUJ 



+ Dk 2 c - E(w) 



(25) 



, R{lu) = v/(l-^i c ) 2 +4a 2 . 

(26) 

The sign in front of the square root R(co) in Eq. (26) agrees 
with the perturbative solution (16). 

To summarize, both the lowest order signal (19) and the 
SCBA are perturbative expressions. However, the lowest or- 
der approximation (19) is valid for a <C 1, while utilizng the 
SCBA allows us to work up until a < 1 (as described be- 
low). This "convergence enhancement", while questionable 
mathematically, appears to be quite useful practically, as we 
demonstrate below both by comparing the self-energy (26) to 
the Monte Carlo simulations, and by applying it to interpret 
proton spectra in human blood. 



Comparison with Monte Carlo simulations 

In Fig. 2 we compare the above results with Monte Carlo 
simulation of diffusion and relaxation in the 2d medium de- 
scribed in Fig. 1. The numerical self-energy was calculated 
according to Eq. (1) after adjusting central frequency of s(uj) 
for a small shift (ft) due to higher-order correlators. 

Practically, the perturbative self-energy (9) agrees perfectly 
with numerics for a < 0.3. For these values of a the char- 
acteristic "triangular" shape of Re E(w) deviates from simple 
Lorentzian (16) due to contribution of harmonics with k < k c . 
Interestingly, for intermediate a ~ 0.5, the shape of E(w) be- 
comes qualitatively closer to the Lorentzian (16). Indeed, for 
larger a, spins dephase before they can explore the scales ex- 
ceeding l c , which increases the relative contribution of large-Zc 
harmonics. 

Quantitatively, the SCBA ansatz (26) is notably better than 
both (9) and (16) for a > 0.5. One can attain a perfect agree- 
ment of SCBA with the data by allowing a and t c to be fitting 
parameters; at a > 1, the SCBA fit values gradually deviate 
from the geniune parameters, as illustrated in Fig. 2(d). This 
is a consequence of the properties of the SCBA signal (1), (26) 
in the complex plane of z = uit c : When a > 1, the simple 
pole of (1) at z = —io? dives under the branch cut connecting 
z = —i±2a. Physically, the developed perturbative approach 
must break down for a > 1, as the system crosses over from 
the diffusion-narrowing to static dephasing regime, a 3> 1. 
The signal in the latter limit coincides with the characteris- 
tic function of the probability distribution of the local Larmor 
frequency, s(t) = /e~* n ( r )*). The connection of s(uj) to the 
mesoscopic structure in this limit was studied only for dilute 
suspensions of objects with basic geometries [20, 21, 32-36]. 



Comparison with experiment 

We now apply our general results to interpret the line shape 
of water protons in blood as measured by Bj0rnerud et al. 
[27]. Blood plasma was titrated with a superparamagnetic 



(a) 



(b) 



0.2 



2r 





FIG. 2: Self-energy E(ui) obtained in Monte Carlo simulations (thick solid lines, red: Re E; blue: Im E), compared with the perturbative result 
for Re E pcrt [Eq. (9) with correlator ^(q) from Fig. 1] (black dashed). Also shown results for the simplified model (13) using approximation 
(16) with properly chosen k c [cf. Fig. 1] (dash-dotted red and blue), the SCBA (26) (dash-dotted magenta and cyan), and the fit of ReE to 
SCBA ansatz (26) (dashed magenta directly on top of Re E), with free parameters t c and a. The values of coupling constant (18) are a = 0.23 
(a), a = 0.79 (b), and a = 1.0 (c). Panel (d) shows the ratio of SCBA fit parameters if\ a fit , and 8Q Rt = a Rt /t^ to their genuine values as 
function of the true a = SO, ■ t c . 



contrast agent to match the magnetic susceptibility of deoxy- 
genated red blood cells. The mesoscopic magnetic structure 
originated from the susceptibility contrast between plasma 
and oxygenated haemoglobin in erythrocytes. 

Our model and the expression (26) for the lineshape has 
been obtained assuming unrestricted diffusion (uniform dif- 
fusivity D). Strictly speaking, the assumption of homoge- 
neous diffusivity does not hold for blood. Below we par- 
tially account for the hindered diffusion in blood by using the 
reduced apparent diffusion coefficient (ADC) of blood [37]. 
This is a reasonable approximation in view of the fast ex- 
change through the cell membrane [38, 39]. 

We have calculated the corresponding self-energy from the 
data of Ref. [27] in three stages: First, the imaginary part 
Ims(w) was determined using the Kramers-Kronig relations 



from the published Res (a;). Special attention was paid to 
phase correction procedure that allowed to cancel the residual 
uniform Larmor frequency offset. Second, the self-energies 
E oxy and Edeoxy for the oxygenated and deoxygenated states 
respectively were determined according to Eq. (1). Finally, 
the difference E(w) = E oxy — Edeoxy was formed in order to 
cancel microscopic effects and reduce data processing errors. 

As demonstrated in Fig. 3, the functional shape (26) agrees 
very well with the self-energy E(cj). The parameters from fit- 
ting ReE(w) are a fit = 1.3, tf = 1.5 ms, yielding <5fi fit w 
0.86 ms _1 , and l c as 2.9 fj,m assuming the blood ADC value 
D = 1.1 fim 2 /ms at the temperature T = 37 C of the mea- 
surement [27] (we calculated this ADC value from the mean 
diffusivity of the blood ADC measured at 25 C [37] asuming 
Arrhenius temeprature dependence similar to that of water). 
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FIG. 3: Fit of the human blood spectrum from Ref. [27] (solid lines, 
red: Re E, blue: Im E) to the SCBA self-energy (26) (dashed). 





Fitting of Im £ gives results within 5%. The value 5fl is in a 
reasonable agreement with the experimental Sfl = 0.69ms _1 . 
The value of l c matches the size of the doughnut-shaped ery- 
throcytes with large diameter 7 /im and thickness about 2 /im. 



DISCUSSION 



FIG. 4: (a) The symbolic representation of the Born series (30). 
Wavy lines represent interaction with the static heterogeneous Lar- 
mor frequency offset — iQ.^. Solid lines represent the free propaga- 
tors GVk (6). (b) The series (31) for the ensemble-averaged propa- 
gator in terms of the self-energy, (c) Perturbative contributions to the 
self-energy. The first diagram is given by Eq. (9). (d) Symbolic form 
of the SCBA Eq. (24). 



In this work we suggested the general representation ( 1 ) for 
spectral lineshape s(u>) in the diffusion-narrowing regime, in 
terms of the self-energy that contains information about 
mesoscopic relaxation. We underscore that, while nominally 
s(uj) is measured, it is the quantity E(w) that characterizes 
the mesoscopic medium, in a sense that it is trivial when the 
medium is magnetically homogeneous. We related the self- 
energy dispersion to the structural characteristics of magnetic 
medium. 

The present treatment clearly illustrates the challenges of 
quantifying magnetic media below the spatial resolution: The 
self-avaraging property of the measurement implies that two 
media are equivalent from the point of the MR signal if their 
correlation functions coincide. As the sensitivity to the higher- 
order correlators T„ drops fast with increasing n, it is the 
lowest order T n , in particular, n = 2, that are most impor- 
tant. Since there are an infinite variety of media whose lowest 
order T n coincide, the inverse problem is, strictly speaking, 
unsolvable. As with any ill-posed inverse problem, having 
prior knowledge about the system is crucial. In particular, 
knowledge about the number of characteristic length scales in 
the geometric profile of susceptibility allows one to select the 
minimal number of the basis functions (14) and to construct 
the corresponding SCBA ansatz by generalizing the form (26). 
Applying the proposed approach in biomedical MRI may al- 
low one to quantify biophysical tissue properties that can be 
further related to physiological processes and malfunctions. 

We conclude by noting that the present approach based on 
the lineshape (1) and on the simple form of the correlator (13) 



is completely general. With straightforward modifications, it 
can be applied to "resolve" the mesoscopic details in diffu- 
sion, in conductivity, and in light propagation in heteroge- 
neous condensed matter systems. 
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APPENDIX: METHODS 

The calculation of the averaged propagator G is done in two 
stages: (i) finding the propagator Q(r, ro;t) of Eq. (3), and (ii) 
averaging over the realizations of 51 (r). Below we describe 
these stages making use of symbolic notation adopted from 
quantum theory [25]. 

On the stage (i), one uses the fact the exact Green's func- 
tion Q = £ _1 is the inverse of the Bloch-Torrey differential 
operator 



C = d t - DV 2 -U, U 



-in. 



(27) 



or, equivalently, CQ 
bare propagator G° — 



= <5(r — ro)5(t). We now define the 
Cn 1 as the fundamental solution of the 
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diffusion equation, 

C G" = 5(r - v )6(t) , C = d t - DV 2 . (28) 

The function G°, whose Fourier transform is Eq. (6), has a 
familiar form in three dimensions, 

G°(r; t) = 9(t){A7rDt)-^ 2 e- r ^ 4Dt (29) 

[where 6(t > 0) = 1 and 9(t < 0) = is the step func- 
tion]. The exact Green's function Q is then obtained by the 
summation of the operator geometric series 

g = (c - uy 1 = (1 - G° * u)- 1 * G° 

= G° + G° *U *G° + G° *U *G° *U *G° .. (30) 

The latter series is analogous to the Born series for the scat- 
tering amplitude in quantum theory [25, 40]. The Born series 
can be schematically represented by the sum of the Feynman 
graphs (Fig. 4). 

The distribution-averaging stage (ii) formally amounts to 
substituting the products of fi(ri)...fl(r„) by the correspond- 
ing correlators T n . We denote this by joining the crosses in 
the Feynman diagrams (Fig. 4) into all possible combinations 
symbolizing T n . 

The exact averaged propagator in Eq. (4) can be represented 
in terms of the self-energy E Wj k which is a sum of all irre- 
ducible contributions to G (by irreducible diagram we mean 
the one which cannot be cut into two by removing any in- 
ternal solid line) [25]. The self-energy E can be then used 
for the block summation [Fig. 4(b)] to obtain the distribution- 
averaged Green's function 

GV + ..., (31) 

which is equivalent to Eq. (5). 

We end this part by outlining the properties of the exact 
Green's function (5). First, we note that the magnetization 
conservation for short times, s(0) = 1, fixes the normalization 
for large frequency behavior irrespectively of the medium: 

Gu, ,kL— >oo = — - — • (32) 

Second, due to causality, as for any response function, 
G(t, r)| t <o = 0, which requires that G w ,k [and thereby s(w)] 
be analytic in the upper-half-plane of the complex variable lu. 
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